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Abstract 



Steady simple shear flow of a low-density binary mixture of inelastic 
smooth hard spheres is studied in the context of the Boltzmann equation. This 
equation is solved by using two different and complementary approaches: a 
Sonine polynomial expansion and the direct simulation Monte Carlo method. 
The dependence of the shear and normal stresses as well as of the steady gran- 
ular temperature on both the dissipation and the parameters of the mixture 
(ratios of masses, concentration, and sizes) is analyzed. In contrast to previ- 
ous studies, the theory predicts and the simulation confirms that the partial 
temperatures of each species are different, even in the weak dissipation limit. 
In addition, the simulation shows that the theory reproduces fairly well the 
values of the shear stress and the phenomenon of normal stress differences. 
On the other hand, here we are mainly interested in analyzing transport in 
the homogeneous shear flow so that, the possible formation of particle clusters 
is ignored in our description. 
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I. INTRODUCTION 



Many features associated with dissipation in rapid granular flows can be well represented 
by a fluid of hard spheres with inelastic collisions. In the simplest model the grains are taken 
to be smooth so that the inelasticity is characterized by means of a constant coefficient of 
normal restitution. The essential difference with respect to normal fluids is the absence of 
energy conservation, which leads to modifications of the usual hydrodynamic equations. In 
recent years, the Boltzmann and Enskog equations have been generalized to account for 
inelastic binary collisions. These equations have been solved by means of an expansion akin 
to the Chapman- Enskog method up to the Navier-Stokes order and detailed expressions for 
the corresponding transport coefficients have been obtained These expressions are 

not restricted to the low-dissipation limit and comparison with Monte Carlo simulations 
indicate that the results are very accurate, even for strong dissipation p. In the context 
of multicomponent granular gases, most of the existing work appears to be based on weak 
dissipation approximations Given that the inelasticity is small, an usual assumption 

in these studies is to consider a single temperature variable characterizing the entire mix- 
ture. However, as one of the authors pointed out P, the equipartition of energy is not 
completely justified beyond the low-dissipation limit and it is necessary to offer theories 
involving mixtures of granular materials in which the kinetic temperatures of species Tj 
are different from the mixture temperature T. As a matter of fact, recent experiments |§] 
and simulations on driven granular mixtures show that the two types of grains do not 
attain the same granular temperature. In terms of the mean square velocities of species, 
this implies a violation of the classical equipartition theorem. A related finding for a binary 
mixture undergoing homogeneous cooling (i.e., an unforced system) has been reported by 



Garzo and Dufty [p!0| , pA| from a kinetic theory analysis. 

All the above works refer to near equilibrium situations. Very little is known about far 
from equilibrium states. This is true for both molecular and granular fluids due to the intri- 
cacy of the Boltzmann and Enskog collision operators. Nevertheless, the difficulties are even 
harder for granular gases since gradients in the system can be controlled by dissipation in 
collisions and not only by the boundary and initial conditions. Thus, for instance, a granular 
system with uniform boundaries at constant temperature develops spatial inhomogeneities 



One of the simplest far from equilibrium physical situations corresponds to the simple 
shear flow. Macroscopically, it is characterized by uniform density and temperature and a 
constant velocity profile. In the case of molecular fluids, this state is not stationary since the 
temperature increases monotonically in time due to viscous heating. However, for granular 
fluids a steady state is possible when the viscous heating is exactly balanced by the inelastic 
cooling. As a consequence, for a given shear rate, the temperature is a function of the 
restitution coefficient in the steady state. This steady state is precisely what we want to 
analyze here. 

In the case of a one-component system, the simple shear flow has been extensively studied. 



Thus, Lun et al. [O] obtained the rheological properties of a dense gas for small inelasticity. 



while Jenkins and Richman |]T4| used a maximum-entropy approximation to solve the Enskog 



equation. An extension of the Jenkins and Richman work to highly inelastic spheres has 
been recently made [|l6|,|l7l. For low-density granular gases, Sela et al. have been able to 
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get a perturbation solution of the Boltzmann equation to third order in the shear rate, finding 
normal stress differences at this level of approximation. On the other hand, some progresses 



have been done by using model kinetic equations in the low-density limit |T8| as well as 
for dense gases [1^. In both works, comparison with Monte Carlo simulations shows an 
excellent agreement even for strong dissipation. Similar studies for multicomponent systems 
are more scarce. Most of them are based on a Navier-Stokes description of the hydrodynamic 
fields [§-0] and, therefore, they are restricted to small shear rates, which for the steady shear 
fiow is equivalent to the low-dissipation limit. As said before, although these studies permit 
different temperatures for the two species, they lead to equal partial granular temperatures 
Tj in the quasielastic limit. A primary attempt to include temperature differences was 
made by Jenkins and Mancini , although applications of this theory which appear in the 



literature incorporate the assumption of equipartition of energy [^]. Since this assumption 
is not completely justified [ll0| , ^2| -|25|| for highly inelastic spheres, the problem of describing 
the simple shear fiow from a multi-temperature theory is still open. 

The aim of this paper is to get the rheological properties of a binary granular mixture 
subjected to the simple shear fiow in the framework of the Boltzmann equation. Two 
complementary routes are followed. First, the set of coupled Boltzmann equations are 
solved by using a first-Sonine polynomial approximation with a Gaussian measure. The 
main characteristic of our solution is that the reference Gaussian distributions are defined in 
terms of the kinetic temperatures Tj instead of the mixture temperature T. Consequently, 
we do not assume a priori the equality of the three temperatures and the temperature ratio 
T1/T2 is consistently determined from the solution to the Boltzmann equations. It is found 
that the partial temperatures of each species are clearly different and so, the energy is not 
equally distributed between both species. The consequences of this effect on the rheological 
properties are significant, as shown below. Once the temperature ratio is known, we get 
explicit expressions for the elements of the pressure tensor. The results are general and 
no limited to weak inelasticity or specific values of the parameters of the mixture. As a 
second alternative and to test the reliability of the theoretical predictions, we have used the 
Direct Simulation Monte Carlo (DSMC) method |^ to numerically solve the Boltzmann 
equation in the simple shear fiow. Although the DSMC method was originally devised for 
molecular fiuids, its extension to dealwith inelastic collisions is easy p7|^ . For the elements 
of the pressure tensor the agreement between theory and simulation turns out to be very 
good over a wide range of values of the restitution coefficients, mass ratios, concentration 
ratios and size ratios. It must be noted that in this paper we are interested in analyzing 
transport properties in the uniform shear fiow. As several authors have shown, the 



simple shear fiow is unstable to long enough wavelength perturbations so that clusters of 
particles are spontaneously developed throughout the system. Here, we will restrict ourselves 
to the uniform case, assuming that the system has reached such a state, and without paying 
attention to the possible formation of particle clusters (microstructure). 

The plan of the paper is as follows. In Sec. |I|, the coupled Boltzmann equations and the 
corresponding hydrodynamic equations are recalled. The steady shear fiow problem is also 
introduced in Sec. while the Sonine approximation is discussed in Sec. |T|. Section 
deals with the Monte Carlo simulation of the Boltzmann equation particularized for steady 
simple shear fiow. The comparison between theory and simulation is presented in Sec. 
and we close the paper in Sec. |V1| with a short discussion. 
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II. THE BOLTZMANN EQUATION AND THE SIMPLE SHEAR FLOW 



We consider a binary mixture of smooth hard spheres of masses mi and m2 and diame- 
ters (Ti and (72. The inelasticity of colhsions are characterized by three independent constant 
coefficients of normal restitution an, 0:22, and ai2 = Q.21, where aij is the restitution co- 
efficient for collisions between particles of species i with j. In the low-density regime, the 
distribution functions fi{r,v,t) {i = 1,2) for the two species verify the set of nonlinear 
Boltzmann equations |T^ 



(5i + Vi-V)/,(r,Vi,t)=^J,,[vi|/,(t),/,(t)] . (1) 

J 

The Boltzmann collision operator Jij [vi|/j, fj] describing the scattering of pairs of particles 



IS 



Jij [Vl|/i, fj] = Crfj J dV2 J dBQ{a- gl2)(CT ■ gi2) 

X «r.2/^(r,v;,t)/,(r,v^,t)-/,(r,Vi,t)/,(r,V2,t)] , (2) 

where a^j = (a, + aj) /2, a is a unit vector along their line of centers, G is the Heaviside 
step function, and gi2 = Vi — V2. The primes on the velocities denote the initial values 
{v'^, V2} that lead to {vi, V2} following a binary collision: 

v'l = vi - (1 + a^r.i) (5- . gi2)5-, V2 = V2 + fiij (1 + a^^^ {a- ■ gi2)o- , (3) 

where fiij = rrii/ {rrii + rrij). The relevant hydrodynamic fields are the number densities rii, 
the flow velocity u, and the temperature T. They are deflned in terms of moments of the 
distributions /j as 



J rfvi/i(vi) , pu = J2piUi = J2j dvimiVifi{vi) , (4) 



nT = E^^^^ = E/ ^vi^V^iV.(vi) , (5) 

i i 

where n = ni + n2 is the total number density, p = pi + P2 = rniUi + m2n2 is the total mass 
density, and Vi = vi — u is the peculiar velocity. Equations and (^) also deflne the flow 
velocity Uj and the partial temperature Tj of species i. 

The collision operators conserve the number of particles of each species and the total 
momentum, but the total energy is not conserved: 

/rfviJ,,[vi|/„/,] = 0, (6) 
J2 J dviTTiiViJijlvilfi, fj] = , (7) 
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^ J rfviim,V^iV,,[vi|/„ /,] = -fnTC , (8) 

where ( is identified as the "coohng rate" due to inelastic colhsions among all species. At a 
kinetic level, it is also convenient to discuss energy transfer in terms of the "cooling rates" 
(i for the partial temperatures Tj. They are defined as 

= E / rfvilm,y,V.,[vi|/„ /,] . (9) 

The total cooling rate ( can be written as 

( = T-'Y.xmCi, (10) 

i 

Xi = Ui/n being the molar fraction of species i. 

From Eqs. (§)-(§), the macroscopic balance equations for the mixture can be obtained. 
They are given by 

DtUi + n^V ■ u + = , (11) 

rrii 



Dtu + p-^V -P = , (12) 

AT - - E ^ + 7^ (V • q + P : Vu) = -C T . (13) 

In the above equations, Df = dt + u ■ V is the material derivative, 

ii = m,J rfviVii^,(Vi) (14) 
is the mass flux for species i relative to the local flow, 

P = E = E / ^^ViVi /,(Vi) (15) 

i i 

is the total pressure tensor, and 

q = E = E Jd^i Imy^^x /.(vo (le) 

i i 

is the total heat flux. The partial contributions to the pressure tensor, Pj, and the heat flux, 
qj, coming from species i can be identified from Eqs. (|15D and ([T6|). 



As said in the Introduction, here we are interested in evaluating the rheological properties 
of a granular binary mixture subjected to the simple shear flow. From a macroscopic point 
of view, this state is characterized by a constant linear velocity profile u = Uj = a ■ r, where 
the elements of the tensor a are akt = a^kx^iy, a being the constant shear rate. In addition, 
the partial densities and the granular temperature T are uniform, while the mass and 
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heat fluxes vanish by symmetry reasons. Thus, the (uniform) pressure tensor is the only 
nonzero flux in the problem. On the other hand, the temporal variation of the granular 
temperature arises from the balance of two opposite effects: viscous heating and dissipation 
in collisions. In the steady state both mechanisms cancel each other and the temperature 



remains constant. In that case, according to the balance energy equation (|T3|), the shear 
stress Pxy and the cooling rate ( are related by 

3 

aPxy = -i^Cp, (17) 

where p = riT is the pressure. Our aim is to analyze this steady state by means of an 
(approximate) analytical method as well as by performing Monte Carlo simulations of the 
Boltzmann equation. 

The simple shear flow becomes spatially uniform when one refers the velocities of the 
particles to a frame moving with the flow velocity u: /j (r, vi) — > /j(Vi). Consequently, the 
corresponding Boltzmann equations (|I]) read 

- aV,,y^h{Y,) = Jn[Vi|A,/i] + M'Vi\fij2] , (18) 

- aV,,y^f2{Vl) = J22[Vi|/2,/2] + J2l[Vi|/2,/l] . (19) 

The elements of the partial pressure tensors Pj (i = 1, 2) can be obtained by multiplying the 
Boltzmann equations by Vi^kVi,£ and integrating over Vi. The result is 



0'kmPl,em + a.lmPl,km — ^ki + ^ki (1^2) , (20) 



where 



A^, = m,J dV^V^,kVl,eJ^A^l\f^, /,] • (21) 
From Eq. (^), in particular, one obtains 

aPi,xy = -fpiCi , (22) 

0'Pl,yy = ^xy + ^xy ) (23) 

^ ^ ^yy ^ + ^zz ■ (24) 

Here, pi = n\T\ = {Pi^xx + Pi,yy + Pi,zz)/'^ is the partial pressure of species 1 and upon 
writing Eq. (^) we have considered the relation (|[). The corresponding equations for P2 
can be easily written just by interchanging the indices 1 and 2. Thus, the determination 
of the elements of the partial pressure tensors Pj is a closed problem once the cooling rates 
Q and the coUisional moments A^^^ are known. This requires the explicit knowledge of the 
velocity distribution functions /j. 
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III. APPROXIMATE SOLUTION 



Unfortunately, solving the Boltzmann equations (|TB|) and (ITW) is a formidable task and 



it does not seem possible to get the exact forms of the distributions /j, even in the one- 
component case. A possible way to overcome such a problem is to expand fi in a complete 
set of polynomials with a Gaussian measure and then truncate the series. In practice, 
Sonine polynomials are used. This approach is similar to the usual moment method for 
solving kinetic equations in the elastic case. In the context of granular gases, this strategy 
has been widely applied in the past few years in the one-component case as well as for 
multicomponent systems and excellent approximations have been obtained by retaining only 
the first two terms. Therefore, one can expect to get a reasonable estimate for Q and A]^^ 
by using the following approximation for fi: 



in ■ / 1 \ 1 

1 + (v^,kVi,e - g^'^fe^jj , (25) 



where fi^M is a Maxwellian distribution at the temperature of the species i, i.e., 

,,„,vo^„,(^)"%.p(-^). ,.e) 

As we will show later, in general the three temperatures T, Ti, and T2 are different in the 
inelastic case. For this reason we choose the parameters in the Maxwellians so that it is 
normalized to Ui and provides the exact second moment of fi. The Maxwellians fi^M for the 
two species can be quite different due to the temperature differences. This aspect is essential 
in our two-temperature theory and has not been taken into account in all previous studies. 
The coefficient Cj can be identified by requiring the moments with respect to Vi^jJ/i/ of the 
trial function (^) to be the same as those for the exact distribution fi. This leads to 

Pi . . 

CiM = — ■ Ske- (27) 

Pi 

With this approximation, the integrals appearing in the expressions of Q and A^j^^ can be 
evaluated and the details are given in Appendices and 0. 

In order to express the solution of the system of equations for the pressure tensor, it is 
convenient to introduce dimensionless quantities. Thus, we introduce the reduced cooling 
rates (* = Q/i', the reduced temperature T* = z/^/a^, and the reduced pressure tensors P* = 
Pi/xip. The (reduced) total pressure tensor P* = P/p = xiP* -|- X2P2- Here, z/ = ^¥71(7^2^0 



is a characteristic collision frequency and Vq = y2T(mi + m2)/mim2 is a thermal velocity 
defined in terms of the temperature of the mixture T. Notice that, for given values of the 
parameters of the mixture, T* and P* are functions of the restitution coefficients Uij only. 

According to the symmetry of the problem, Pi^xz = Pi,yz = 0, so that the nonzero elements 
are Pi^xx, Pi,yy, Pi,zz, and Pi,xy = Pi,yx- The three normal elements are not independent since 
Pi,xx + Pi,yy + Pi,zz = 37j, where the temperature ratios ji = Ti/T are given by 

7 1 
^1 = ^ 1 TT' ^2 = — J TT , (28) 
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with 7 = T1/T2. The temperature ratio 7 provides information about how the kinetic 
energy is distributed between both species. In addition, according to Eq. (^41) , P*yy = P*^^- 
Consequently, the partial pressure tensors have four relevant elements, say for instance: 
V = {Pi^yy, P2,yyj Pi,xyy P2,xy}- Taking into account the results derived in the Appendices, 
Eqs. ([23| ) and ( pi]) (plus their corresponding counterparts for species 2) can be written as 



CV=Q, 



(29) 



where C is the 4x4 matrix 



1 


-{Gil + G12] 
—H21 




1 

—H12 
-{G22 + G21, 



-(Gil + Gi2)z//a 
-H2iiy/a 





-Hi2u/a \^ 
-{G22 + G2i)iy/a 


/ 



(30) 



and 



Q 



I 








-^11 + P12 
\ F22 + F21 / 



(31) 



Here, the functions Fj,, Gj,, and ifj, are defined in the Appendix 0. The solution to Eq. 



IS 



(32) 



Equation ( p2[ ) gives the nonzero elements of the pressure tensors P.* in terms of the reduced 
temperature T* (or the reduced shear rate a/v), the temperature ratio 7, the restitution 
coefficients and the parameters of the mixture. The dependence of T* on the coefficients 
can be obtained from the energy balance equation (HT 



-1/2 



2 P* 



3 xi7iCr + 3^272(2 
2 xiPt^y + X2-P. 



(33) 



2,xy 



Finally, when Eqs. (|32|) and (^) are used in Eq. ( ^2]) (or its counterpart for the species 2), 
one gets a closed equation for the temperature ratio 7, that can be solved numerically. In 
reduced units, this equation can be written as 



7 



C2 Pi 



xy 



QP 



(34) 



2,xy 



In the elastic limit 



a 



1), we recover previous results derived for molecular gases 
A simple and interesting case corresponds to the case of mechanically equivalent particles 
(mi = 1712, «ii = «22 = «i2 = = (^22)- In this limit, Eqs. (|32D -(^4D leads to 7 = 1, 



p* 



P^, with 



p* 
yy 



22 + a 
3 3 — a 



(35) 
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(37) 



and 

r- . 4 = 5 (i±i!)P^(i _ . (38) 

5 2 + a ^ ^ ^ ^ 



These expressions differ from tlie results derived in Ref. |JT8|| by using a model kinetic equa- 



tion. However, for practical purposes, the discrepancies between both approximations are 
quite small, even for moderate values of the restitution coefficient. It is also interesting to 
consider the limit of weak dissipation (1 — aij <^ 1), in which case it is possible to get ana- 
lytical results. For the sake of simplicity, let us assume that all the particles have the same 
coefficient of restitution, namely, an = = ai2 = To get the first order corrections in 
the quasielastic limit, we introduce the perturbation parameter e = (1 — a^)^/^ and perform 
a series expansion around e = 0. The details of the calculation are presented in Appendix 
and here we only quote the final results. First, the leading term of the reduced shear rate 
a/u (which is a measure of the steady granular temperature) is 

ao{l-aY^ + ---. (39) 



V 

Next, the temperature ratio and the relevant elements of the (partial) pressure tensors can 
be written as 

7 = l + 7oa2(l-«2) + ..., (40) 
i^:,, = l + P/|),a2(l-a2)^...^ (41) 

^% = l + ^S«o(l -«') + ■■■, (42) 

^:.. = ^S«o(l-«^)^/^ + ---. (43) 

In these equations, ao, 7o, and P/^^ are dimensionless coefficients that depend on the ratios 
of mass, concentrations and sizes. Their explicit expressions are given in Appendix y. 

In summary, by using the Sonine approximation (^), we have explicitly determined the 
rheological properties of the mixture as well as the reduced shear rate and the temperature 
ratio as functions of dissipation and mechanical parameters of the mixture. These theoretical 
predictions will be compared with those obtained from Monte Carlo simulations in Sec. |V[ 
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IV. MONTE CARLO SIMULATION 



From a numerical point of view, the Direct Simulation Monte Carlo (DSMC) method ||2^ 



is the most convenient algorithm to study non-equilibrium phenomena in the low-density 
regime. It was devised to mimic the dynamics involved in the Boltzmann collision term. The 
extension of the DSMC method to deal with inelastic collisions is straightforward ||T9| , ^ , ^ 



and here we have used it to numerically solve the Boltzmann equation in the simple shear 
flow. In addition, since the simple shear flow is spatially homogeneous in the local La- 
grangian frame, the simulation method becomes especially easy to implement. This is an 
important advantage with respect to molecular dynamics simulations. On the other hand, 
the restriction to this homogeneous state prevents us from analyzing the possible instability 
of simple shear flow or the formation of clusters or microstructures. 

The DSMC method as applied to the simple shear flow is as follows. The velocity 
distribution function of the species i is represented by the peculiar velocities {V^} of iVj 
"simulated" particles: 

f,(y,t)-.n,-J2S{V-Vk{t)). (44) 

k=i 

Note that the number of particles Ni must be taken according to the relation N1/N2 = ni/n2. 
At the initial state, one assigns velocities to the particles drawn from the Maxwell-Boltzmann 
probability distribution: 

/,(V, 0) = n-'/' K.t3(0) exp (-V'XM , (45) 



where VJ]^(0) = 2T(0)/mj and T(0) is the initial temperature. To enforce a vanishing initial 
total momentum, the velocity of every particle is subsequently subtracted by the amount 
iV~^ J2k Vfc(O). In the DSMC method, the free motion and the collisions are uncoupled over 
a time step At which is small compared with the mean free time and the inverse shear rate. 
In the local Lagrangian frame, particles of each species (z = 1, 2) are subjected to the action 
of a non-conservative inertial force Fj = —rrii a • V. This force is represented by the terms 
on the left-hand side of Eqs. (|IBp and (0). Thus, the free motion stage consists of making 
Vfc — i> Vfc — a ■ VfcAt. In the collision stage, binary interactions between particles of species 
i and j must be considered. To simulate the collisions between particles of species i with 
j a sample of ^NiCU^J^At pairs is chosen at random with equiprobability. Here, uj^2l 
an upper bound estimate of the probability that a particle of the species i collides with a 
particle of the species j. Let us consider a pair {k,i} belonging to this sample. Here, k 
denotes a particle of species i and i a particle of species j. For each pair {k, i} with velocities 
{Vfc, V^}, the following steps are taken: (1) a given direction a-^i is chosen at random with 
equiprobability; (2) the collision between particles k and i is accepted with a probability 

equal to 6(gH ■ 5-H)a;i7V'^miL where u^P = 4:'^<^ijnj\gke ■ o-^l and gM = - V^; (3) if 
the collision is accepted, postcollisional velocities are assigned to both particles according 
to the scattering rules: 

Vfc Vfc - fiji{l + aij){gke ■ o-ki)^ki , (46) 
^ + + aij){gke ■ S-ke)a-k£ . (47) 



10 



In the case that in one of the colhsions to^f > cj^i^, the estimate of w^i^ is updated as 
'^max = ^kP- The procedure described above is performed for i = 1,2 and j = 1, 2. 

In the course of the simulations, one evaluates the total pressure tensor and the partial 
temperatures. They are given as 

P = E^i:V.V,, (48) 

i=i k=i 



^^ = ^EVL (49) 

'^^^i k=l 

To improve the statistics, the results are averaged over a number A/" of independent real- 
izations or replicas. In our simulations we have typically taken a total number of particles 
N = N1 + N2 = 10^ a number of rephcas Af = 10, and a time step At = 3 x 10"3Aii/1/oi(0), 
where An = (y^nniaf^)^^ is the mean free path for collisions 1-1. 

A complete presentation of the results is complex since there are many parameters in- 
volved: {aij, mi/ 1712, ni/n2, (Til/ a22}- For the sake of concreteness, henceforth we will con- 
sider the case an = 022 = cn2 = «• In the steady state, the reduced quantities T*, 7^, 
and P* are independent of the initial state for given values of the restitution coefficient and 
the ratios of mass, concentration and sizes. To illustrate it, in Fig. |l| we present the time 
evolution of T*(t) for a = 0.75, 0"ii/(j22 = 1, mi/m2 = 4, ni/n2 = 1/3 and three different 
initial conditions. Time is measured in units of Aii/Voi(0). After an initial transient pe- 
riod, all curves converge to the same steady value, as predicted by the solution described in 
Sec. fT|. The same qualitative behavior has been found for the temperature ratio and the 
elements of the reduced pressure tensor. Therefore, in the following we will focus on the 
dependence of the steady values of the reduced quantities on the restitution coefficient a 
and the parameters of the mixture, once we have checked they do not depend on the initial 
state. 



V. COMPARISON BETWEEN THEORY AND MONTE CARLO SIMULATIONS 

In this Section we compare the predictions of the Sonine approximation with the results 
obtained from the DSMC method. Our goal is to explore the dependence of a*, 7 = T1/T2 
and the nonzero elements of P* on a, the mass ratio fi = mi/m2, the concentration ratio 
5 = ni/n2, and the ratio of sizes w = CTii/a22- 

First, we will investigate the dependence of the relevant quantities on a and fi for given 
values of 6 and w. Recent molecular dynamics simulations for a dilute monocomponent 
system of smooth inelastic hard disks |31| have supported an "equation of state" to a sheared 
granular system in which the steady (reduced) temperature T* can be closely fitted by a 
linear function of (1 — Similar results have been obtained from kinetic models of the 



Boltzmann and Enskog |jT9| equations. An interesting question is whether this simple 
relationship can be extended to the case of multicomponent systems. The results obtained 
here (both from the simulations and from the kinetic theory analysis) for mixtures of different 
masses, concentrations or sizes show that T* is indeed a quasi-linear function of (1 — a^)~^. 
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As an illustrative example, we consider the case w = 1, S = 1 (equimolar mixture), and three 
different values of the mass ratio /i = 1, 2, and 10. Figure |^ shows T* versus (1 — as 
obtained from the simulations (symbols) and from the Sonine approximation (lines). It is 
evident that the kinetic theory has an excellent agreement with the simulation results and 
also that T* is practically linear in (1 — a^)^^ whatever the mass ratio considered is. The 
slope of the straight lines increases as the disparity of masses increases. 

The temperature ratio and the non-zero elements of the pressure tensor are plotted in 
Figs. H and ^ (a — ci), respectively, as a function of the dissipation parameter a for the same 
cases as those considered in Fig. ^ The curves corresponding to // < 1 can be easily inferred 
from them. Figure ^ clearly shows that, except for mechanically equivalent particles, the 
partial temperatures are different, even for moderate dissipation (say a ^ 0.9). This means 
that the traditional assumption of equipartition of fluctuation energy begins to fail. This 
effect is generic of multicomponent dissipative systems and is consistent with results recently 
derived in the homogeneous cooling state |TD|,^ as well as in driven systems |]^,|^ . The extent 
of the equipartition violation depends on the concentrations and the mechanical differences of 
the particles (e.g., masses, sizes, restitution coefficients), and is greater when the differences 
are large. The agreement between theory and simulation is again excellent. With respect 
to the pressure tensor. Fig. | (a — ci), we see that the theory captures well the main trends 
observed for the rheological properties. At a quantitative level, the agreement is better in 
the case of the shear stress P*^ and the normal element P*^, while the discrepancies for the 
normal stresses P*y and P*^ are larger than in the case of the temperature ratio, especially 
as the restitution coefficient decreases. On the other hand, the theory only predicts normal 
stress differences in the plane of shear flow {P*^ ^ P*y = P*^) while the simulation also 
shows that there is anisotropy in the plane perpendicular to the flow velocity, P*^ > P*y. 
This kind of anisotropy has been also observed in molecular dynamics simulations of shear 
flows ||3^. Nevertheless, these relative normal stress differences in this plane are very small 
and decrease as a increases. 

The influence of the concentration ratio 6 on the temperature ratio and the rheological 
properties is shown in Figs. |^ and ^ (a — rf), respectively, for w = 1, /i = 4, and two values 
of 6: 5 = 1/3 and 6 = 3. We observe again a strong dependence of the temperature ratio 7 
on dissipation. For a given value of a, the temperature ratio increases as the molar fraction 
of the heavy species decreases. Concerning shear stresses, we see that they are practically 
independent of the concentration ratio since all the curves collapse in a common curve. A 
more signiflcant influence is observed for the normal stresses. In general, the agreement with 
the theory is good although the discrepancies are more important in the case of 5 = 1/3. 
Finally, the influence of the size of the particles on the rheological properties is illustrated in 
Fig. ^ (a — (i). We consider an equimolar mixture [6 = 1) of particles of equal mass (/i = 1) 
for two different values of the size ratio: w = 1, and w = 2. We see that similar conclusions 
to those previously found in Figs. 1^-^ are obtained when one considers mixtures of particles 
of different sizes. 



VI. DISCUSSION 

In this paper we have addressed the problem of a low-density granular mixture consti- 
tuted by smooth inelastic hard spheres and subjected to a linear shear flow = ay. We 
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are interested in the steady state where the effect of viscosity is compensated for by the 
dissipation in colhsions. Our description apphes for arbitrary values of the shear rate a or 
the inelasticity of the system and not restriction on the values of masses, concentrations 
and sizes are imposed in the system. The study has been made by using two different and 
complementary routes. On the one hand, the set of coupled Boltzmann equations are solved 
by means a Sonine polynomial approximation and, on the other hand, Monte Carlo simula- 
tions are performed to numerically solve the Boltzmann equations. Given that the partial 
temperatures Tj of each species can be different, the reference Maxwellians in the Sonine 
expansion are defined at the temperature for that species. This is one of the new features of 
our expansion. On the other hand, to put this work in a proper context, it must be noticed 
that we have restricted our considerations to states in which the only gradient is the one 
associated with the shear rate so that density and velocity fluctuations are not allowed in 
the numerical simulation. 

We have focused on the analysis of the dependence of the steady (reduced) temperature 
T* and the (reduced) pressure tensor P* on the coefficients of restitution aij and the pa- 
rameters of the mixture, namely the mass ratio /i, the concentration ratio 5 and the size 
ratio w. The results clearly indicate that the deviation of the above quantities from their 
functional forms for elastic collisions is quite important, even for moderate dissipation. In 
particular, the temperature ratio, which measures the distribution of kinetic energy between 
both species, is different from unity and presents a complex dependence on the parameters 
of the problem. This result contrasts with previous results derived for granular mixtures 
where it was consistently assumed the equahty of the partial temperatures in the small 



inelasticity limit. In the same way as in the homogeneous cooling state problem |IO| , p8 
the deviations from the energy equipartition can be weak or strong depending on the me- 
chanical differences between the species and the degree of dissipation. On the other hand, 
the simulation as well as the theoretical results also show that the steady total temperature 
T* can be fitted by a linear function of (1 — a^)~^ with independence of the values of the 
parameters of the mixture. This extends previous results derived in the context of simple 



granular gases by using molecular dynamics ||3T| or Monte Carlo simulations [|T8ip!9i . With 
respect to the rheological properties, comparison between theory and simulation shows a 
good quantitative agreement, especially for the shear stress P*y, which is the most relevant 
element of the pressure tensor in a shearing problem. Although the kinetic theory also pre- 
dicts normal stresses, the discrepancies between theory and simulation are larger than those 
found for the temperature ratio or the shear stress. 

It is illustrative to make some comparison between the predictions made from our two- 
temperature theory with those obtained if the differences in the partial temperatures were 
neglected (Ti = T2 = T). For instance, let us consider the mixture an = (722, ni = n2, 
and mi = 10m2 with a = 0.75. In this case, for the xy and yy elements of the pressure 
tensor, the simulation results are —P^y = 0.498 and P*y = 0.723. Our two-temperature 
theory predicts —P^y = 0.498 and P*y = 0.743 while the single-temperature theory (as- 
sumption made in previous works) gives —P^y = 0.456 and P*y = 0.815. Clearly, inclusion 
of the two-temperature effects improves the theoretical estimations and makes a significant 
(quantitative) difference with respect to the predictions of the single-temperature theory. 

As a final comment, let us mention that the study made here can in principle be extended 
in both aspects, kinetic theory and simulations, to the revised Enskog equation in order to 
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assess the influence of finite density on the rheological properties of the mixture. Work along 
this line is in progress. 
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APPENDIX A: EVALUATION OF THE COOLING RATES 



In this Appendix the (reduced) cooling rates (* are evaluated by using the first Sonine 
approximation (^). The cooling rate is given by 



(Al) 



where 9i = = Vi/vq, J*- = {vl/ninal2)Jij, and /* = {vl/ni)fi. Henceforth, 

it will be understood that dimensionless quantities will be used and the asterisks will be 
deleted to simplify the notation. A useful identity for an arbitrary function /i(Vi) is given 
by 



dVi/i(Vi)J.,[Vi|/„/,] =^o{^^ J dVi J f/V2/.(Vi)/,(V2) 



(A2) 



with 



V'l = Vi - + aij){B ■ gi2)o- . 



(A3) 



This result applies for both i = j and i ^ j- Use of this identity in Eq. CKl) allows the 
angular integrals to be performed. The result is 







:i - 
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-V 



(^12/ 



rfVi / rfV2(/?2/^(Vl)/.(V2) 



+ {l-al)^^e,iilx^ j dVi j rfV2^?2/.(Vi)/,(V2) 

+ {l + a,,)'^^9,^l,,x, j rfVi j rfV2(7i2(gi2-G,,)/,(Vi)/,(V2), (A4) 

where Gij = fJ.ij'Vi + fJ,jiV2- Now we consider the Sonine approximation (^) for the distri- 
butions fi'. 



/.(Vi 



71 



3/2 



1 + 



(A5) 



Neglecting nonlinear terms in the tensor Ci^ke^ the expression (|A4|) can be written as 
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dVi / rfV2(?i2(gi2-G,,) e-(^' 



(A6) 



Here, use has been made of the fact the scalar Q cannot be coupled to the traceless tensor 
Ci^ke so that the only nonzero contribution to the cooling rate comes from the Maxwellian 
term (first term of the right hand side of (|A5|) ) of the distribution function. The first integral 
of Eq. ( |A6| ) is straightforward and can be done with a change of variables to relative and 
center of mass variables. The next two integrals are somewhat more complicated and they 
can be performed by the change of variables 



Vi - V2, y = ^^.Vi + 9,\2, 



(A7) 



with the Jacobian {6i + 6j) ^. The integrals can be now easily performed and the final result 
for ^1 is 



Ci 



0"12/ 



"11 y 



'91 + 02 



1/2 



(l + ai2) 



2 - /i2i(l + ai2. 



91 + 02 
02 



(A8) 



The result for ^2 is obtained from Eq. 



by interchanging 1 and 2. 



APPENDIX B: EVALUATION OF THE COLLISIONAL MOMENTS 



In reduced units, the collisional moments A^"'^ are given by 



rriiV, 



T 

rUiVQ 



^-1/2 



dYiVi,kVi,iJ.j[Vi\f„fj 



T ^ 



1/2 f ^ij \ ' 



Vcri2/ 

x(a- ■ gi2) (VikVii - Vi^kVi,, 



rfVi / rfV2/,(Vi)/,(V2) 



where the identity ( |A2| ) has been used. Substitution of (^) into Eq. 
integral to be performed with the result 



da-e{d- ■ gi2) 

(Bl) 

allows the angular 



+ / t/Vi / rfV2/,(Vi)/,(V2) 



X 



2 T Vai2/ 

gi2{Gij,kgi2,i + Gij^eQu^k) + ^(3 - aij)gi2gi2,k9i2/ - ^(1 + aij)gl2^, 



(B2) 



where (712,^ = ^i,fc — ^2,^ and Gij^k = fi'ijyi,k + fi'jiV2,k- To perform the integral we use the 
Sonine approximation of and the change of variables (|A7| ) . When one neglects again 



nonlinear terms in the tensor Cj, the collisional moment A^j^^ becomes 
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X 

X 



(712/ ' {Oi + ej) 

1 + OiiOi + ej)-^Q ■ (y + ^,x)(y + ^,x) + 9^(9^ + O^yKj ■ (y - ^ix)(y - ^,x) 



{6i + 6j) ^x{xkye + xeUk) + XijXXkXe — ^(1 + aij)x^Ske 



(B3) 



where 



bij = 6i6j{6i + 6j) ^, 



(B4) 



dij — {9i + 9j) ^, 



(B5) 



(B6) 



The corresponding integrals can be now easily performed and, after some algebra, the final 
result is 



3/2 



X 



A 9 ^^^^ 
Okl — I- 



1 + lA 



(B7) 



Prom this general expression one can get the coUisional moments Al\, Alj, A^^^, and A^J. 
In particular. 



(B8) 



(B9) 



where 



i^ii = :i^V^/^2"M (^)%r'^'(l + an)(2 + aii), 

io VCl2/ 



(BIO) 



V (Ji2 / 



(Bll) 



3/2 



5A12 + 2A*2l(l + 0:12) 



(B12) 
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1/2 



'3^ 9^ + 62 



(B13) 



1/2 



do 



\ 3^ ^^1 + ^2^ 



(B14) 



The corresponding expressions for F22, 6*22, -^21, G21, and H21 can be easily inferred from 
Eqs. (|ET0| )-(|ET4l) by just making the changes 1^2. From Eqs. (|B|)-(|BT|) and (^), it is 
easy to prove the identity 



7iCi = -^11 + F12 + (Gil + G12) 7i + -f^i272, 



(B15) 



which is in fact required by the partial energy conservation equation (|2^) to support the 
solution found for the simple shear flow problem. In the case of mechanically equivalent 
particles, the expression of the coUisional moments A^^^ are 



A'^, + Al] = —{l + a 



'■kl 
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[2 + a)6M--iS-a)PM 



(B16) 



When a = 1, Eq. ( |B16| ) reduces to the results derived from the Boltzmann equation in the 
first Sonine approximation. 



APPENDIX C: LOW DISSIPATION LIMIT 



In this Appendix we derive the expressions for the main quantities of the simple shear 
flow problem in the low-dissipation limit. The expansions in powers of e = (1 — a^)^/^ are 
given by Eqs. (|39|)-(^3D. For symmetry reasons, the expansion of P*^y has only odd powers, 
while those of the normal stresses (and of the temperature ratio) have only even powers. 
However, from a practical point of view, it is simpler to use a* = a/u as a perturbation 
parameter instead of e. Thus, the expansions (|39| ) -(|43| ) can be written in terms of a* as 



a 



1 + ctoa*^ + 



(CI) 



7 = 1 + 7oa*^ + 



P* =14- P.(^) 0*"^ -\ 



P* = 1 + P.^^^ fl*^ 4- ■ ■ 

l,XX ' l,XX^ ' 



(C2) 
(C3) 
(C4) 



p* 

i,xy 



(C5) 



Of course, both expansions are directly related, so that 
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ao 



V-2ao' 



(C6) 



The coefficients Pj^^/}^ can be obtained from Eq. (|32D by retaining terms up to second 
order in a*. To do that, we need the corresponding expansions of the quantities Fij, Gij, 
and Hij. They are given by 



Fn 



^11 + -Til a + ■ 



(2)^*2 



;|ao - |a;27o)a*^ H 



(C7) 



Gil = -|y2^a:i (— (l + '^-foX2a*' + ■ ■ ■) 



5 ' ■ \o-i2/ 



16 / 5ao + 70(9x2 - 7/i: 

-/i2lX2 (^1 + 



F -^^ 
^12 — — 



■12) „*2 



(C8) 



(C9) 



G 



12 



16 
15 



2;2yU2i(3 + 2/i 



, 5ao/^i2 - 7o [4/i?2 + Aii2(2a;i - 1) + 3(a;i - 1)] ^2 , 



G<°> 



'^12 



+ 



(CIO) 



H 



12 



16 

15 



X2l^l2 



-/i2i(ao - 270/^12)0* + ■ ■ 



-^12 + -"12 « + 



(Cll) 



The expressions of the coefficients F^p, G[^\ and H-j^ can be easily identified from the above 
equations. On the other hand, the coefficients F22, G22, F21, G21, and if2i are obtained by 
just making the changes 1 ^ 2, ao ^ «o and 70 —70- Substitution of Eqs. (|C7|) -( |CTID 
(and their counterparts for species 2) into Eq. ( ^2] ) and considering only terms through 
second order in a* allows one to get the coefficients P/^^^ in terms of ao and 70. The final 
expressions of such coefficients will be omitted here since they are very large and not very 
illuminating. 

Once the pressure tensors are known, we are in conditions to get the coefficients ao and 
7o. First, the cooling rates behave as 



cr 



cTii Y 



'2fi2i [^j ^1 + 2yU2ia;2 j ao - 4/i2iAii2a^27o 
(Aiao + 5i7o) a*2 + . . . , 



+ . . 



(C12) 



'2/ii2 y—j ^2 + 2Aii2Xi I ao + 4/ii2/i2ia;i7o 



{A2ao + B2I0) a*' + -- 



(C13) 
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where again the coefficients Ai and Bi are easily identified. The quantities ao and 70 can be 
now easily obtained from the partial balance of energy (p2|), i.e., 



-lp,% = Aao + Bao (^ = 1,2). (CM) 



3 

The solution of this system of algebraic equations leads to 



"° = 3 AA-AA ■ 



where 



p(i) 



_ 2 A2Pi ly — AlP2jy 



(C16) 



Hu {F22 



(0) 



^(0) , ^(0)^ 
<-^12 + '-^21 . 



.(0) 



.(0) 



(o)^2 I 77(0)77(0)1 

+ -"12-^21 



'^21 



(Gil + <-^12 Jl<^22 + <^21 ) - tiyi ti; 



.(O)vr^(o) 



(o)^ 



r(0) rr(0) 
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(C17) 



and P2^y can be obtained from P\\xy by setting the changes 1^2. 



(1) 



and the coefficients Pfli are obtained from Eq. ([ 
in the dissipation parameter e are retained. 



In summary, Eqs. ( |C15D and ( |C16| ) give ao and 70, respectively, while is given by ( |(J6D 

when only terms through second order 
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FIGURES 



FIG. 1. Time evolution of the reduced granular temperature T*{t) = iy^{t)/a^ as obtained 
from Monte Carlo simulation of the Boltzmann equation for a — 0.75. an = a22, ni/n2 = 1/3, 
mi/m2 = 4, and starting from three different initial conditions. Time is measured in units of 
Aii/yoi(0). 

FIG. 2. Plot of the reduced granular temperature T* = ja? versus the parameter (1 — ct^)^^ 
as obtained from simulation (symbols) and the Sonine approximation (lines), for w = o\\l(j22 = 1; 
6 = ni/n2 = 1 and three different values of the mass ratio = ■mi/m2- /x = 10 (solid line and 
circles); n = 2 (dashed line and squares), and /x = 1 (dotted line and triangles). 

FIG. 3. Plot of the temperature ratio 7 = T1/T2 as a function of the restitution coefficient a 
as obtained from simulation (symbols) and the Sonine approximation (lines). We have considered 
"W = o'ii/(T22 = 1, 5 = ni/n2 = 1 and three different values of the mass ratio /i = mi/ni2: fi = 10 
(solid line and circles); fi = 2 (dashed line and squares), and = 1 (dotted line and triangles). 

FIG. 4. Plot of the reduced elements of the pressure tensor (a) —Pxy = —Pxy/p-, (b) 
Pxx — Pxx/Pi (c) Pyy = Pyy / P and (d) P*^ = Pzz/p versus the restitution coefficient a as 
obtained from simulation (symbols) and the Sonine approximation (lines). We have considered 
w = 0'ii/c22 = 1, S = ni/n2 = 1 and three different values of the mass ratio /i = mi/ 1112: /i = 10 
(solid line and circles); fj, = 2 (dashed line and squares), and = 1 (dotted line and triangles). 

FIG. 5. Plot of the temperature ratio 7 = T1/T2 as a function of the restitution coefficient a 
as obtained from simulation (symbols) and the Sonine approximation (lines). We have considered 
w = crii/o'22 = 1, M = mi/m2 = 4 and two different values of the concentration ratio S = ni/n2- 
(5 = 1/3 (solid line and circles), and 6 = 3 (dashed line and squares). 

FIG. 6. Plot of the reduced elements of the pressure tensor (a) —Pxy = —Pxy/p, (b) 

P*x = Pxx/p, (c) Pyy = Pyy/p and (d) P*^ = Pzz/p versus the restitution coefficient a as 
obtained from simulation (symbols) and the Sonine approximation (lines). We have considered 
w = o\\/a22 = 1, M = mi/m2 = 4 and two different values of the concentration ratio 5 = ni/n2'- 
6=1/3 (solid line and circles), and 6 = 3 (dashed line and squares). 

FIG. 7. Plot of the reduced elements of the pressure tensor (a) —Pxy = —Pxy/p, (b) 
Pxx = Pxx/p, (c) Pyy = Pyy/p and (d) P*^ = Pzz/p versus the restitution coefficient a as 
obtained from simulation (symbols) and the Sonine approximation (lines). We have considered 
fj, = nil/ 171,2 = 1, 6 = ni/n2 = 1, and two different values of the size ratio w = Cw/a22'- w = 2 
(solid line and circles), and w = 1 (dashed line and triangles). 
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